Ultrafast effective multi-level atom method for primordial hydrogen recombination 
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Cosmological hydrogen recombination has recently been the subject of renewed attention because 
of its importance for predicting the power spectrum of cosmic microwave background anisotropies. 
It has become clear that it is necessary to account for a large number n > 100 of energy shells 
of the hydrogen atom, separately following the angular momentum substates in order to obtain 
sufficiently accurate recombination histories. However, the multi-level atom codes that follow the 
populations of all these levels are computationally expensive, limiting recent analyses to only a 
few points in parameter space. In this paper, we present a new method for solving the multi-level 
atom recombination problem, which splits the problem into a computationally expensive atomic 
physics component that is independent of the cosmology, and an ultrafast cosmological evolution 
component. The atomic physics component follows the network of bound-bound and bound-free 
transitions among excited states and computes the resulting effective transition rates for the small 
set of "interface" states radiatively connected to the ground state. The cosmological evolution 
component only follows the populations of the interface states. By pre-tabulating the effective 
rates, we can reduce the recurring cost of mufti- level atom calculations by more than 5 orders of 
magnitude. The resulting code is fast enough for inclusion in Markov Chain Monte Carlo parameter 
estimation algorithms. It does not yet include the radiative transfer or high-n two-photon processes 
considered in some recent papers. Further work on analytic treatments for these effects will be 
required in order to produce a recombination code usable for Planck data analysis. 

PACS numbers: 98.80.Es, 98.62.Ra, 32.80.Ec 



I. INTRODUCTION 

The advent of high precision cosmic microwave back- 
ground (CMB) experiments, such as Planck has re- 
cently motivated several authors to revisit the theory of 
cosmological recombination pioneered by Peebles [2j and 
Zeldovich et al. Q in the 1960s. The free electron frac- 
tion as a function of redshift x e (z) is one of the major 
theoretical uncertainties in the prediction of the CMB 
temperature and polarization anisotropy power spectra 
To obtain a recombination history accurate to 
the percent level, it is necessary to account for a high 
number of excited states of hydrogen, up to a princi- 
pal quantum number n max = 0(100) 0,1a]. The desired 
sub-percent accuracy can only be reached when explicitly 
resolving the out-of-equilibrium angular momentum sub- 
states, which requires the multi-level atom (MLA) codes 
to follow A^icvd = n max (n max + l)/2 individual states. 
Moreover, the ordinary differential equations (ODEs) de- 
scribing the level populations are stiff, requiring the so- 
lution of large Mevci x N\ cvc \ systems of equations at each 
integration time step. This problem has been solved by 
several authors [9hLU| , but each of these codes takes hours 
to days to run. 

Eventually, it is necessary to be able to produce not 
only accurate but also fast recombination histories, to be 
included in Markov Chain Monte Carlo (MCMC) codes 
for cosmological parameter estimation. The MCMC re- 
quires CMB power spectra (and hence recombination his- 
tories) to be generated at each proposed point in cosmo- 
logical parameter space, with a typical chain sampling 
0(1O 5 ) points iH]. Furthermore, dozens of MCMCs 
are often run with different combinations of observa- 



tional constraints and different parameter spaces. This 
makes it impractical to include recombination codes 
that run for more than a few seconds in the MCMC. 
One solution is to precompute recombination histories 
x e (z\H , T CM b, £l m h 2 ,tt b h 2 , Yho, N„) on a grid of cosmo- 
logical parameters, and then use elaborate interpolation 
algorithms to evaluate the recombination history for any 
cosmology [lH, or to construct fitting functions pi Il4{. 
However, such procedures need to be re-trained every 
time additional parameters are added, and are rather un- 
satisfying regarding their physical significance. 

In this work we present a new method of solution for 
the recombination problem, perfectly equivalent to the 
standard MLA method, but much more efficient compu- 
tationally. The basic idea is that the vast majority of 
the excited hydrogen levels are populated and depopu- 
lated only by optically thin radiative transitions (bound- 
bound and bound-free) in a bath of thermal photons; we 
show that their effect can be "integrated out" leaving 
only a few functions of the matter and radiation temper- 
atures T m and T T (this list would include the free elec- 
tron density n e if we incorporated collisions), which can 
be pretabulated. In an actual call to the recombination 
code from an MCMC, it is then only necessary to solve an 
effective MLA (hereafter, EMLA) with a smaller number 
of levels (perhaps only 3: 2s, 2p and 3p), which elimi- 
nates the computationally difficult Aicvci x Ai cve i system 
solution in the traditional MLA. [The idea is similar in 
spirit to the line-of-sight integral method for the com- 
putation of the CMB power spectrum , which elimi- 
nated a large number of independent variables from the 
cosmological perturbation theory system of ODEs (the 
high-order moments of the radiation field, Of for £ ^> 1) 
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in favor of pretabulated spherical Bessel functions.] Our 
method achieves a speed-up of the recombination calcu- 
lation by 5 to 6 orders of magnitude. 

We note that our method only eliminates the compu- 
tational complexity associated with the high-n excited 
states and does not include continuum processes and ra- 
diative transfer effects that have been studied by pre- 
vious authors p^| - [26j . However, we note that there has 
been much progress in analytic treatments of these effects 
[21M23I ] ; ultimately, we expect to improve these analytic 
treatments and graft them (and an analytic treatment of 
helium recombination |27H3l| ) on to the ultrafast EMLA 
code described herein to yield a recombination code that 
is accurate enough for Planck data analysis. 

This paper is organized as follows. In Section [TT] we 
review the general picture of hydrogen recombination, 
and the bound-bound and bound-free transition rates 
involved in the calculation. In Section IIIII we describe 
the standard MLA method. We present our new EMLA 
method in Section IIVI and demonstrate its equivalence 
with the standard MLA formulation. We describe our 
numerical implementation and results in Section |Vl and 
conclude in Section I VII Appendix [S] is dedicated to 
demonstrating the invertibility of the system defining the 
EMLA equations. Appendix IBl proves a complementar- 
ity relation between effective transition probabilities. We 
prove detailed balance relations between effective transi- 
tion rates in Appendix [Cl Appendix [Dl exposes the post- 
saha approximation we use at early times when comput- 
ing recombination histories. 



II. BOUND-BOUND AND BOUND-FREE 
TRANSITION RATES 

The evolution of the free electron fraction is governed 
by the network of transitions between bound states of 
hydrogen as well as recombination and photoionization 
rates. Before giving detailed expressions for these rates, 
let us first outline the general picture of the process of 
recombination. 

It has long been known that direct recombinations to 
the ground state are ineffective for recombination 0, Q , 
since the resulting emitted photons are immediately reab- 
sorbed by hydrogen in the ground state, as soon as the 
neutral fraction is higher than ~ 10~ 9 . Electrons and 
protons can efficiently combine only to form hydrogen 
in excited states. The minute amount of excited hydro- 
gen at all relevant times during cosmological recombina- 
tion is not sufficient to distort the blackbody radiation 
field near the ionization thresholds of the excited states. 
Recombination to the excited states is therefore a ther- 
mal process: it depend on the matter temperature T m 
which characterizes the free electrons and protons veloc- 
ity distribution, and also on the radiation temperature 
T r , since the abundant low-energy thermal photons can 
cause stimulated recombinations. Photoionization rates 
from excited states only depend on the radiation tem- 



perature since they do not involve free electrons in the 
initial state. 

Transitions between bound excited states may be ra- 
diative or collisional. Radiative transition rates are well 
known and depend only on the radiation temperature 
characterizing the blackbody radiation field, undistorted 
in the vicinity of the optically thin lines from the Balmer 
series and beyond. Collisional transition rates are much 
less precisely known, but only depend on the matter tem- 
perature and the abundance of charged particles caus- 
ing the transitions (i.e. free electrons and free protons, 
which, once helium has recombined, have the same abun- 
dance n e — ?ip due to charge neutrality). 

Finally, some of the excited states can radiatively de- 
cay to the ground state. The most obvious route to 
the ground state is through Lyman transitions from the 
p states. However, due to the very high optical depth 
of these transitions, emitted Lyman photons are imme- 
diately reabsorbed by hydrogen atoms in their ground 
state. This "bottleneck" can only be bypassed by the 
systematic redshifting of photons, which can escape re- 
absorption once their frequency is far enough below the 
resonant frequency of the line. The relevant transition 
rate in this case is the net decay rate to the ground 
state, which is a statistical average accounting for the 
very small escape probability of Lyman photons. Two- 
photon transitions are usually much slower than single- 
photon transitions. However, the rate of two-photon de- 
cays from the metastable 2s state is comparable to the 
net decay rate in the highly self-absorbed Lyman transi- 
tions, and this process should therefore be included in a 
recombination calculation 0, Q • 

We now give explicit expressions for the bound-bound 
and bound-free rates dicussed above. Subscripts nl re- 
fer to the bound state of principal quantum number n 
and angular momentum quantum number I. We denote 
af s the fine structure constant, /i e = m e m p / (m e + m p ) 
the reduced mass of the electron-proton system, Ej the 
ionization energy of hydrogen, and E n = —Ejn^ 2 the en- 
ergy of the n th shell. Finally, we denote by fBB{E,T r ) = 
(e £ / Tl — I) -1 the photon occupation number at energy 
E in the blackbody radiation field at temperature T r . 



A. Recombination to and photoionization from the 
excited states 

The recombination coefficient to the excited state 
nl, including stimulated recombinations, is denoted 
c*ni(T m , T r ) (it has units of cm 3 s _1 ). The photoioniza- 
tion rate per atom in the state nl is denoted /? ra ;(T r ). 
Both can be expressed in terms of the bound-free radial 
matrix elements g(n,l, k,1') (32j |. Defining 

lM = ^c4^(l + nV) 3 

x J2 mz*(l,l')g(n,l,K,l') 2 , (I) 
i'=i±i 
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where k denotes the momentum of the outgoing electron 
in units of h/ag (where do is the reduced-mass Bohr ra- 
dius), the recombination coefficient is given by (32[: 



a„i (T m , T T 



h 3 



(27r Me T m ) 3 /2 

" + OO 



E ' K/T -lm{n) 
x{l + f BB (E Kn ,T T )}d(K 2 ), (2) 



x / e 
lo 



where E Kn = Ej(k 2 + n~ 2 ). The photoionization rate 
only depends on the radiation temperature and can be 
obtained by detailed balance considerations from the re- 
combination coefficient: 



, N (27TU e T r ) 3 / 2 p /T 



Tr,T T ) 



(3) 



B. Transitions between excited states 

We denote R n i^n'i' the transition rate from the excited 
state nl to the excited state n'V , It has units of sec -1 
per atom in the initial state. Transitions among excited 
states can be either radiative or collisional: 

Rnl^n'V — R T nl^>n>l' C^r) + ^n°-m'/' C^mi n e)i (4) 

where n e = n p is the abundance of free electrons or free 
protons. In this paper, we follow exclusively the radiative 
rates. These are given by 



R 



rad 

nl—tn'l' 



l' [1 + f&B{E nn < , T T j] E n > E n 



2ltp- E n>n/ T t fi» ad 
gi c "nT-mi 



E n < E n i 



(5) 



where E nn ' = E n — E n > is the energy difference between 
the excited levels, gi = 21 + 1 is the degeneracy of the 
state nl, and A n i, n 'i' is the Einstein ^-coefficient for the 
nl — > n'V transition, which may be obtained from the 
radial matrix element i?"L|33l: 



2tt 3 Ei 



1 



I 



max(Z, I') 
21 + 1 



r>nl 

■tU.ni 



I |2 



(G) 



C. Transitions to the ground state 



In the Sobolev approximation [34] with optical depth 
T np.is 3> 1; the net decay rate in the np —> Is transi- 
tion is: 



«np\ lg 



8ttH 



ip •^ x lsfnp) 
{ x np ~ Sxisfnp) , (7) 



where X n = he/ E n \ is the transition wavelength, and 
/+, is the photon occupation number at the blue side of 
the corresponding Ly-n line. In this paper, we will take 
fn P = fB&{E n i,T r ), i.e. assume the incoming radiation 
on the blue side of the line has a blackbody spectrum. 
This assumption is actually violated due to feedback from 
higher-frequency Lyman lines (e. g. radiation escaping 
from Ly/3 can redshift into Lya) [27l l35l . |36| ; while our 
formalism is general enough to incorporate different 
we have not yet implemented this in our code. 

The 2s state cannot decay to the ground state through 
a radiatively allowed transition. This decay is how- 
ever possible with a two-photon emission, which, al- 
though slow, is comparable in efficiency to the highly 
self-absorbed Lyman transitions. The simplest expres- 
sion for the net 2s — > Is two-photon decay rate is: 



"*2s| ls = A 2s i s \X 2s 



xi s e 



-E 2 /T t 



), (8) 



where A2 S i s ~ 8.22 s" 1 is the total 2s — > Is two-photon 
decay rate [37| . 

In each case, we denote the net downward rate in the 
i — > Is transition, where i 6 {2s, 2p, 3p, ...}: 



xn 



XiRi- 



x ls Rn 



(9) 



where the rates R depend on atomic physics, T r , and the 
optical depths in the Lyman lines. 

Both the Sobolev approximation for the np — > Is tran- 
sitions Eq. ([7]), and the simple expression Eq. ([5]) for 
the net 2s — > Is two-photon decay do not account for 
subtle yet important radiative transfer effects. An accu- 
rate recombination calculation should account for time- 
dependent effects in Lyq l2ll [2^| a suite of two-photon 
continuum process es \yA Il8t fell fe~5j . and resonant scat- 
tering in Lya (23l . |26j. These are not included in the 
present code and we plan to add them in the future us- 
ing analytic treatments. 



Finally, the ground state population a;i s w 1 — x e 
evolves due to transitions from and into the np and 2 s 
states (two-photon transitions from higher energy states 
are dominated by "1 + 1" photon decays, already ac- 
counted for). Photons emitted in the Lyman lines are 
very likely to be immediately reabsorbed, and the only 
meaningful quantity for these transitions is the net de- 
cay rate in the line, which is a statistical average over 
a large number of atoms, and accounts for the very 
low escape probability of a photon emitted in the line. 



III. THE STANDARD MLA METHOD 

Although the standard MLA formulation does not 
make this distinction, we cast the excited states of hy- 
drogen into two categories. On the one hand, most ex- 
cited states are not directly radiatively connected to the 
ground state. We call these states "interior" states and 
denote Xk the fractional abundance of hydrogen in the 
interior state K € {3s, 3d, 4s, 4d, 4/, 5s, ...}. On the other 
hand, the 2s and np states (n > 2) are directly radiatively 
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connected with the ground state. We call these states "in- 
terface" states and denote the fractional abundance of 
hydrogen in the interface state i G {2s, 2p, 3p, ...}. 

In the standard MLA formulation, the free electron 
fraction x e (z) is evolved by solving the hierarchy of cou- 
pled differential equations: for the interior states, 

X K = x 2 e n H a K + ^ x lRl^k + XjRj^K 

L 3 

-x K (p K + J2 r k^l+12 Rk ^)' ( 10 ) 

L 3 

for the interface states, 
it = xlnncti + ^ XlRl^i + 2J x i R j •■ + xuRu-yi 

L j 

-x l (^/3 l + Y,R l ^L + Y, R i^3 + Ri^u); (ii) 

L 3 

and for the free electrons and ground state, 

X e = -its = Xis^Ris^i - y^XjRj^u. (12) 

i i 

The radiative rates between excited states are many or- 
ders of magnitude larger than the rate at which recombi- 
nation proceeds, which is of the order of the Hubble rate. 
Even the relatively small net rates out of the interface 
states (A2s,i s and A2p t i a /T2 Pt u) are still more than 12 or- 
ders of magnitude larger than the Hubble rate. The pop- 
ulations of the excited states can therefore be obtained 
to high accuracy in the steady-state approximation (this 
approximation is ubiquitous in many problems and has 
long been used in the context of cosmological recombi- 
nation @, [1(1 HI|, where its accuracy has been tested 
explicitly [HI). Setting Xk and ii to zero in Eqs. (fTO]) 
and pT]) , we see that the problem amounts to first solve a 
system of linear algebraic equations for the Xk , Xi , with 
an inhomogeneous term depending on x e , and then use 
the populations Xi in Eq. (TT2"j) to evolve the free elec- 
tron fraction. The solution of the system of equations 
(|10p. (jlip needs to be done at every time step, since the 
inhomogeneous term of the equation depends on the ion- 
ization history, which explicitly depends on time as well 
as on the cosmological parameters. Recent work [l(| [H| 
has shown that to compute sufficiently accurate recombi- 
nation histories, one needs to account for excited states 
up to a principal quantum number n max ~ 100, resolving 
the angular momentum substates. This requires solving 
an O (l0 4 x 10 4 ) system of equations at each time step, 
which, even with modern computers, is extremely time 
consuming. 

IV. NEW METHOD OF SOLUTION: THE 
EFFECTIVE MULTI-LEVEL ATOM 

We now give a computationally efficient method of so- 
lution for the primordial recombination problem. We fac- 



tor the effect of the numerous transitions involving inte- 
rior states in terms of effective transitions into and out 
of the much smaller number of interface states. Once the 
rates of these effective transitions are tabulated, the cos- 
mological evolution of the free electron fraction can be 
obtained from a simple effective few-level atom calcula- 
tion. We describe the method in Section IIV Al and give 
the proof of its exact equivalence to the standard MLA 
method in Section llVBl In the subsequent Section ITV CI 
we consider which states should be treated as interface 
states. 



A. Motivations and general formulation 

We first note that the only quantity of importance for 
CMB power spectrum calculations is the free electron 
fraction as a function of redshift, x e (z). The populations 
of the excited states are calculated only as an interme- 
diate step - if they are desired (e.g. to calculate Ha 
scattering features [38]), the populations of the excited 
states can be obtained by solving Eqs. (|T0l ITT]) once the 
free electron fraction is known. Furthermore, only the 
"interface" states 2s and np are directly connected to 
the ground state and directly appear in the evolution 
equation for the free electron fraction Eq. ([T2l . All other 
( "interior" ) excited states are only connected with other 
excited states or with the continuum, through optically 
thin radiative transitions (and to a lesser extent through 
collisions §)■ Interior states are only transitional states: 
an electron in the "interior" rapidly transitions through 
spontaneous and stimulated decays or absorptions caused 
by the blackbody radiation field (or collisions with free 
electrons and protons), until it is either photoionized, or 
reaches an interface state. There can be a very large num- 
ber of transitions before any of these outcomes occurs, 
but the passage through the "interior" is always very 
short compared to the overall recombination timescale, 
and can be considered as instantaneous (for the same 
reason that the steady-state approximation is valid in 
the standard MLA formulation). 

Instead of computing the fraction of hydrogen in each 
interior state K, one can rather evaluate the probabil- 
ities that an atom initially in the interior state K ul- 
timately reaches one of the interface states or gets pho- 
toionized. Of course, after reaching an interface state, the 
atom may perfectly transition back to an interior state, 
or get photoionized. However, we consider the proba- 
bility of first reaching a given interface state before any 
other one, which is uniquely defined. For an atom in 
the interior state K, we denote by P K the probability 
of ultimately reaching the interface state i, and P K the 
probability of ultimately being photoionized. The proba- 
bilities P l K must self-consistently account for both direct 
transitions K — > i and all possible indirect transitions 
K — > L — > i (with an arbitrary number of intermediate 
states). Mathematically, this translates to the system of 
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FIG. 1. Schematic representation of the formulation of the 
recombination problem adopted in this work. Dotted arrows 
represent possibly numerous fast transitions within the "inte- 
rior" . 



linear equations: 



pk = J2^ p i 



R 



K^i 



r 



K 



r 



K 



(13) 



where T k is the total width (or inverse lifetime) of the 
state K: 



K ■ 



(14) 



Similarly, the P% must satisfy the self-consistency rela- 
tions: 



pe 



R 



K- 



L 



r 



-p? 



K 



T k 



(15) 



We show in Appendix [5] that these linear systems are 
invertible and therefore uniquely determine P l K and P^. 
In Appendix |B] we prove the complementarity relation, 



(16) 



which has the simple physical interpretation that an atom 
in the Kt\i interior state eventually reaches an interface 
state or is photoionized with unit probability. 

Once these probabilities are known, it is possible to 
describe the large number of transitions between all the 
states in a simplified manner, in terms of effective rates 
into and out of the interface states. To clarify the expla- 
nation, we illustrate in Figure Q] the processes described 
below. 

An electron and a proton can effectively recombine to 
the interface state i either through a direct recombina- 
tion (with coefficient ccj), or following a recombination 
to an interior state K (with coefficient olk), from which 
a sequence of interior transitions may ultimately lead to 
the interface state i with probability P K . The effective 
recombination coefficient to the interface state i is there- 
fore: 



Ai = ai + ot K P K - 



(17) 



Conversely, an atom in the interface state i may effec- 
tively be ionized either through a direct photoionization 
(with rate ft), or after being first excited to an interior 
state K (with rate R^k), from which the atom may ul- 
timately be photoionized after a series of interior transi- 
tions with probability P K . The effective photoionization 
rate from the interface state i is therefore: 



(18) 



K 



Finally, atoms can effectively transition from an interface 
state i to another interface state j, either through a di- 
rect transition if it is allowed, or after first transitioning 
through the interior. The effective transfer rate between 
the ith and jth interface states is therefore: 



n 



t->3 



Ri 



-^1 



K 



U + i). (19) 



A" 



The rate of change of the population of the interface state 
i is therefore: 

&i = xlnnAi + r,R, ., + xis&is^ (20) 

-XifBi+J^^j+Ri^ls), 

where we have included the effective transitions described 
above, as well as transitions from and to the ground state. 

The system of equations (|13ll20|) is exactly equivalent 
to the standard MLA formulation, as we shall show in 
Section HVBl below. 

Let us now consider the dependences of the effective 
rates. In the purely radiative case, the probabilities 
P K and P K only depend on the radiation temperature 
T r , since transitions between excited states and pho- 
toionizations only depend on the locally thermal radia- 
tion field. As a consequence, the effective recombination 
rates Ai (T m ,T r ) are only functions of matter and radia- 
tion temperatures and the effective photoionization and 
bound-bound rates Bi (T r ) and IZ^j (T r ) are functions 
of the radiation temperature only. When including col- 
lisional transitions, all effective rates become functions 
of the three variables T T ,T m and n e . In all cases, effec- 
tive rates can be easily tabulated and interpolated when 
needed for a recombination calculation. 

Intuitively, we would expect that Ai, Bi, and 7£j_>j 
satisfy the detailed balance relations, 



9i> 



(21) 



and 



K 



gie-W'Bifr) = ^i^U(T m = T T ,T t ). (22) 

We show in Appendix [C] that these equations are indeed 
valid. This means that we only need to tabulate half of 
the 7£j_j.j [the other half can be obtained from Eq. (|21l) ] 
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and all the Ai [the Bi can be obtained from Eq. (|22j) : in 
particular, we do not need to solve for the Pjj-]. 

We note that the probabilities P % k,Pk are a general- 
ization of the cascade matrix technique introduced by 
Seaton [39[. Seaton's calculation assumed a vanishing 
ambient radiation field, so that electrons can only "cas- 
cade down" to lower energy states. In the context of the 
recombination of the primeval plasma, one cannot ignore 
the strong thermal radiation field, and electrons rather 
"cascade up and down," following spontaneous and stim- 
ulated decays or photon absorption events. The spirit 
of our method is however identical to Seaton's cascade- 
capture equations (39|, where the "cascading" process is 
decoupled from the particular process populating the ex- 
cited states, or from the depopulation of the interface 
states. 



B. Equivalence with the standard MLA method 

This section is dedicated to proving the equivalence of 
the EMLA equations, Eqs. P^H2U1) . with the standard 
MLA equations, Eqs. (fTO! ITT]) . in the steady-state limit 
for the interior states (i.e. where we set Xk ~ 0). The 
steady-state approximation does not need to be made 
for the interface states to demonstrate the equivalence 
of the two formulations (but we do use it for practical 
computations since it is valid to very high accuracy). 

We denote by A the number of interior states and n* 
the number of interface states (we will address in Section 
IIV CI the issue of which states need to be considered as 
interface states). 

We begin by defining the A x A rate matrix M whose 
elements are 



M KL = 5 KL T K - (1 - 5 KL ) R k ^l- 



(23) 



We also define the n* + 1 length- A vectors P J , P e whose 
elements are the probabilities P % K and P^ respectively, 
and the n* + 1 length- A vectors R 4 , R e of components 

R-k = R-K-^i, (24) 
R e K = 0k. (25) 

The defining equations for the probabilities, Eqs. (1131 
IT51) . can be written in matrix form MP* = R* and 
MP e = R e respectively (after multiplication by IV). 
We show in Appendix lAlthat the matrix M(T r ) is invert- 
ible, for any temperature T r > 0. The formal solutions 
for the probabilities are therefore 



P' 
P< 



M _1 R* 
M _1 R e 



(26) 
(27) 



We also define the length-A vector X which contains the 
populations of the interior states Xk, and the length-A 
vector S of components 



A careful look at Eq. (fTOl) in the steady-state approxi- 
mation (Xx — 0) shows that it is the matrix equation 
M T X = S, which has the solution: 

X = (M T )~ 1 S = (M- 1 ) T S. (29) 
Both Eqs. ([TTj) and (|2T))) can be cast in the form 



->i + x lsRls~¥i 
Ri—tls J Xi I interior- 



(30) 



The only a priori different term is the net transition rate 
from the interior to the state i, ii\ interior- In the standard 
MLA formulation, Eq. (|11|) . this term is 



iL^rtl = E ( X * R K ~ X ^K) (31) 
K 

= X T R 4 

K 



Xi 2J Ri^K- 



(32) 



With our new formulation, Eq. (|20p , using the definitions 
of the effective rates Eqs. pTrfl9)) . the net transition rate 
from the interior to the state i is: 



, (EMLA) 



E 



K 



x e n H a K Pk 
x l R l ^ K {P'k 



x jRj^KPk 



E^) 



(33) 



Using the complementarity relation Eq. (|16[) . we rewrite 
P k + P K = 1 ~ P k- We then recognize that the 

common factor of P % K is just the if-th component of the 
vector S, Eq. (l28l) . so we can rewrite Eq. (|33t as 

(34) 



I (EMLA) _ Q T 



% I interior 



s ip, 



Xi E Ri^K- 
K 



From the formal solution for the populations of the inte- 
rior states Eq. (|29|) . we have 

(35) 



X T R* = STM-^ 1 = S T P J 



where the second equality is obtained from the formal 
solution for the probabilities P % K , Eq. (f2l)]). We therefore 
see from Eqs. and ([34]) that 



, (MLA) 



. I (EMLA) 
* I interior ' 



(36) 



and hence the two formulations are exactly equivalent. 
They only differ by the order in which the bilinear prod- 
uct S" 1 ]^^ is evaluated. 



Choice of interface states 



s 



K 



xlnnaK + E x jRj 



(28) If one naively includes all np states up to n — n max = 

0(100) in the list of interface states, the interpolation of 
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effective rates can become somewhat cumbersome as it 
invoives (10 4 ) functions of one to three variables. How- 
ever, only the lowest few of these states actually have sig- 
nificant transition rates to the ground state; indeed, most 
of the decays to the ground state proceed through either 
2s (two-photon decay) or 2p (Lya escape), as anticipated 
in the earliest studies 0, Q. 

The rate of Lyman line escape is dominated by the 
lowest few lines. For example, if the relative populations 
of the excited states were given by the Boltzmann ratios 
(which is a good approximation until late times) then the 
net decay rate in the np —> Is transition (not accounting 
for feedback from the next line) would be proportional 
to 



oc (1 



2 ) 3 e _jE " /Tr . 



(37) 



This relation would imply that the Ly/3 escape rate is 
< 1% of the Lya escape rate, and the higher-order Ly- 
man lines contribute even less. Our previous compu- 
tations of the escape rates (e.g. Ref. 21|) agree with 
this expectation. These considerations imply that for 
n > 3, Xu\np <C ii s |2p in Eq. (fT2j). Moreover, an atom 
in the np state with n > 3 is much more likely to spon- 
taneously decay to n's or n'd, with 2 < n! < n, than 
to emit a Lyman-n photon that successfully escapes the 



line. This implies that 



in Eq. (fill . 



•Enp | Is 

In addition to a very low net decay rate out of the np 
states for n > 3, feedback between neighboring lines fur- 
ther suppresses their efficiency as interface states. The 
few photons that escape the Ly(n + 1) line will be re- 
absorbed almost certainly in the next lower line, after a 
redshift interval 



A = z em - z ab = (1 + z a b) ( - 1 



(38) 



Feedback between the lowest-lying lines is not instan- 
taneous: Az/(1 + z a b) = 0.185 for Ly/3 — >Lya feedback, 
0.055 for Ly7 — s-Ly/3, and 0.024 for Ly6 — >\sy^. However, 
for higher-order lines, feedback rapidly becomes nearly 
instantaneous as Az/(1 + z a b) ~ 2/n 3 . Thus the effect 
of the higher Lyman lines is even weaker than Eq. (|37[) 
would suggest. Recent work 36] has shown that includ- 
ing lines above Ly/3 results in a fractional error |Ax e |/x e 
of at most ~ 3 x 10~ 4 . 

We therefore conclude that very accurate recombi- 
nation histories can be obtained by only including 
2s, 2p, ...,n*p as interface states and negelecting higher- 
order Lyman transitions altogether. We will use n* = 3 
in this paper, and investigate the optimal value of this 
cutoff more quantitively in future work. 

Our formulation in terms of effective transition rates 
and interface states is therefore much better adapted for 
a fast recombination calculation that the standard MLA 
formulation. To compute accurate recombination histo- 
ries, explicitly accounting for high-n shells of hydrogen, 
one first needs to tabulate the {Ai} and {TZ^j} on tem- 
perature grids. The computation of the effective rates 



is the time-consuming part of the calculation; however, 
since they are independent of the cosmological parame- 
ters, this can be done once, and not repeated for each 
cosmology. The free electron fraction can then be com- 
puted very quickly for any given cosmology by solving the 
n* + 1 equations (|2"0"|) and (fT2)l . interpolating the effective 
rates from the precomputed tables. Note that Eq. (f2"0")) 
is a simple n* x n* system of linear algebraic equations 
in the steady-state approximation. 



V. IMPLEMENTATION AND RESULTS 

Here we give some details on the implementation of our 
EMLA code. Section IV AI describes the computation of 
the effective rates (the computationally expensive part 
of the calculation, which needs to be done only once). 
Section IV Bl descibes the implementation of the ultrafast 
effective few-level atom calculation. We show our recom- 
bination histories and compare our results with the ex- 
isting standard MLA code RecSparse 10] in Section 

Ec] 



A. Computation of the effective rates 

We have implemented the calculation of the effective 
rates in the purely radiative case. Bound-free rates were 
computed by numerically integrating Eq. ^ using an 
11-point Newton-Cotes method, where the radial matrix 
elements g(n,l, k,V) were obtained using the recursion 
relation given by Burgess (32j. Einstein ^-coefficients 
were computed by using the recursion relations obtained 
by Hey [4fJ| for the radial matrix elements Fi- 
nally, we obtained the probabilities PI- using a sparse 
matrix technique similar to that of Ref. [10j] when solving 
Eq. (flUl) . We accounted explicitly for all excited states up 
to a principal quantum number n max , resolving angular 
momentum substates. We tabulated the effective rates 
Ai(T m ,T T ) on a grid of 200 log-spaced points in T r from 
0.04 to 0.5 eV and 20 linearly spaced points in T m /T r 



from 0.8 to 1.0, and 7£j_^(T r ) on the grid of points in 
T r . The maximum relative change in the effective rates 
(Aln^4i or A In Hi-tj ) over the whole range of temper- 
atures considered is 0.051 when comparing n max = 64 
vs. 128, 0.015 when comparing n max = 128 vs. 250, and 
0.005 when comparing n max = 250 vs. 500. 

In the left panel of Figure [H we show the total effective 
recombination coefficient AB(T m , T T ) = A2 S (T m ,T I ) + 
A.2 P (T m ,T r ) computed for n* = 2 (i.e. with inter- 
faces states 2s and 2p only, neglecting all Lyman transi- 
tions above Lyman-a), normalized to the case-B recom- 
bination coefficient aB(T m ). Note that aB(T m ) is just 
As(T m ,T r — 0|7i max = oc) with our notation; indeed, 
for T r = 0, p K = and therefore P C K = for all K 

so Y.i p k - 1 and nence Hi A = Hi a i + J2k a K = 
Hni a nh where the last sum is over all excited states. 
We can see that as the radiation temperature increases 
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(i.e. as the redshift increases in Figure [2]), the conver- 
gence with n max becomes faster. This is to be expected, 
since for higher T r , highly excited hydrogen is more eas- 
ily photoionized, i.e. P^ d becomes closer to unity. In 
that case adding more shells to the calculation does not 
matter so much because recombinations to the highest 
shells are very inefficient, due to the high probability of 
a subsequent photoionization. 

In the right panel of Figure [2] we show the ratio 
•A-2s{T m , T r ) / AB(T m , T T ), which is the fraction of recom- 
binations to the n — 2 shell that are to the 2s level. This 
fraction is in general different from the intuitive value of 
1/4, and its exact value depends on temperature. 



B. Ultrafast EMLA code 

In order to actually compute the recombination his- 
tory, we require an evolution equation for the free elec- 
tron fraction, 



Xe{x e ,nn,H, T m ,T r ), 



(39) 



and in some cases a similar equation for T m . For con- 
creteness, we implement the case of 3 interface states 
i e {2s,2p,3p} (n* = 3). 

To compute x e , we first obtain the downward TZ^j (T T ) 
from our table via cubic polynomial (4-point) inter- 
polation and Ai(T m ,T T ) via bicubic interpolation (2- 
dimensional in lnT r and T m /T T using 4x4 points). The 
upward lZj^i(T r ) are obtained using Eq. ([21]) and the 
effective photoionization rates Bj(T r ) are obtained using 
Eq. We then solve for the {xi} using Eq. (|20[) . and 

finally obtain x e using Eq. (fi"2|) . 

The matter temperature is determined by the Comp- 
ton evolution equation, 



T m = —2HT n 



8<7Ta T T^x e (T r - T m ) 
3(1 + /ho + x e )m e c 



(40) 



where ctt is the Thomson cross section, a T is the radia- 
tion constant, /h c is the He:H ratio by number of nuclei, 
m e is the electron mass, and c is the speed of light. At 
high redshift, one may use the steady-state solution (see 
Appendix A of Ref. [Hi), 



T, Ri T 



1 



3(1 + /hc + x e )m e cH 
8aTa T T^x e 



(41) 



At the highest redshifts, the ODE describing hydrogen 
recombination is stiff; therefore for z > 1570 we follow the 
recombination history using perturbation theory around 
the Saha approximation, as described in Appendix [D] 
At 500 < z < 1570 we use Eq. (gTJ to set the matter 
temperature, and a fourth-order Runge-Kutta integra- 
tion algorithm (RK4) to follow the single ODE for x e (z); 
and at z < 500 we use RK4 to follow the two ODEs for 
x e (z) and T m (z) simultaneously. The integration step 
size is Az = — 1.0 (negative since we go from high to low 
redshifts). 



C. Results and code comparison 

We have tabulated the effective rates for n max = 16, 
32, 64, 128, 250 and 500. It is in principle possible to 
compute the effective rates for an arbitrarily high n max , 
but it is not meaningful to do so as long as collisional 
transitions are not properly accounted for. The recur- 
ring computation time of our ultrafast EMLA code is 
0.08 seconds on a MacBook laptop computer with a 2.1 
GHz processor, independently of n max . Our recombina- 
tion histories are shown in Figure |31 We compared our 
results with the existing standard MLA code RecSparse 
for n max = 16, 32, 64, 128 and 250. As can be seen in 
Figure |H the two codes agree to better than 8 x 10~ 5 
across the range 200 < z < 1600, despite having different 
methods for accounting for the excited states, and inde- 
pendent implementations for matrix elements and ODE 
integration. 



VI. CONCLUSIONS AND FUTURE 
DIRECTIONS 

We have shown that the computation of primordial hy- 
drogen recombination can be factored into two indepen- 
dent calculations. On the one hand, most excited states 
are not directly radiatively connected to the ground state, 
and undergo transitions caused by the thermal bath of 
blackbody photons at the relevant frequencies, as well as 
the thermal electrons and protons. One can account for 
these numerous transitions with effective transition rates 
into and out of the "interface" states which are connected 
to the ground state. The computationally intensive as- 
pect of a recombination calculation in fact resides in the 
evaluation of these effective rates, which are functions of 
matter and radiation temperature only. This calculation 
being independent of cosmological parameters, it can be 
done prior to any recombination calculation, once and for 
all. A simple effective few-level atom can then be evolved 
for any set of cosmological parameters, without any need 
for "fudge factors" or approximations. 

This work does not present a final recombination code 
satisfying the accuracy requirements for future CMB ex- 
periments. Firstly, collisional transitions were not in- 
cluded. They may be particularly important for the high- 
n states. The effective rates computed here are therefore 
only approximating the correct rates in the limit of zero 
density. Our formalism is general and collisions can be 
included as soon as accurate rates are available (the main 
change would be that the interpolation tables would re- 
quire In n e as an additional independent variable) . Sec- 
ondly, we have not included important radiative transfer 
effects, such as feedback between low- lying Lyman lines 
[35l l3o | , two- photon decays from n > 3 [Til ll9l - [2lT |25| , 
resonant scattering in Lya [22|, [23|, [26| , or overlap of the 
high- lying Lyman lines (work in preparation). To pre- 
serve the computational efficiency of our method, fast 
analytic approximations have to be developed to include 
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FIG. 2. Left panel: "Exact fudge factor" as a function of redshift Ab{T^, T r )/aB(Tm), for several values of n max , using T m (z) 
computed by RecSparse for cosmological parameters as in Ref. We use the fit of Ref. [4l[ for the case-B recombination 
coefficient ctB(T m ). For comparison, the code Recfast uses a constant fudge factor F = 1.14 to mimick the effect of high-n 
shells. Right panel: Fraction of the effective recombinations to the n = 2 shell that lead to atomic hydrogen in the 2s state. In 
both cases the effective rates were computed for n» = 2, i.e. with interface states 2s and 2p only, neglecting escape from the 
Lyman /3,7, ... lines. 




1600 




200 400 



600 800 1000 1200 1400 1600 



FIG. 3. Left panel: Relative differences between recombination histories computed with successively more accurate effective 
rates. Right panel: Recombination history for effective rates computed with n max = 500, i.e. accounting explicitly for 125,250 
states of the hydrogen atom. 
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these effects, which will be the subject of future work. 
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FIG. 4. A comparison of our ultrafast code to RecSparse 
[l0| . for different values of n max . The vertical axis is the 
fractional difference in free electron abundance rescaled by 
10 5 (positive indicating that RecSparse gives a larger x e ). 
We see that the maximum fractional deviation is < 8 x 10~ 5 . 
The feature around z = 1540 is due to a timestep change in 
RecSparse. 
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Appendix A: Invertibility of the system defining the 

Pk,PK- 

In this section we show that the matrix M(T r ) defined 
in Eq. (|23[) is non-singular, for any value of the radiation 
temperature T T > 0. 

Let us consider the eigenvalue equation Mb = and 
select a particular K\ such that I&r-J > for all L. 
The eigenvalue equation implies 

= \M KlKl b Kl + M KlL b L \ 

|- \MK lL \\bL\, (Al) 

where we have used the inverse triangle inequality. The 
matrix M is diagonally dominant, i.e. 

VK, M KK = T k > J2 Rk^l = \ Mkl \ ( A2 ) 

Using the inequality (fMt for K = K\ in Eq. (fATj) . wc 
obtain 



o> l^ild^J-N 



(A3) 



For any interior state , there always exists a sequence 
of transitions that ultimately leads to some interface state 
i, K\ — > K 2 — > ... — > K n — > i, for any temperature 
T r > 0. (i.e. there are no "dead end" interior states). 
In particular, (Mr-^! = Rk 1 ^k 2 > 0. For Eq. (|A3|) 
to hold, it is therefore necessary that fr^J = 1 6 j<- a | . Re- 
peating the above reasoning recursively leads to \bK t \ = 
l^ 2 | = -.. = |^J. 

For the last interior state of this sequence, K n , the 
inequality (|A2[) is strict since Rn n ^i > 0. The eigenvalue 
equation projected on K n leads to Eq. (|A1|) for K n : 



(A4) 



> M KnKn \b Kn \ - \ M K n L\\b L \- 

If b ^ 0, then \bx n \ > and the strict inequality (|A2[) 
for K — K n used in Eq. (|A4[) implies the contradictory 
result 



0> \ M k~l\ 



> 0. 



(A5) 



As a consequence, Mb = implies that b = neces- 
sarily. This proves that M(T r ) is nonsingular, for any 
T r > 0. 



Appendix B: Proof of the complementarity relation 

Pk + Pk = ^ 

We define the length- N vector V = (1, 1, 1) T , and 
note that 

(MV)jc - ]T Mkl = Yl Rk ^ + P« ( B1 ) 



(the Rk^l terms cancel). In matrix form, this reads: 



MV = Y RJ 



R e 



M 



(B2) 



The matrix M being invertiblc, this implies J^j P J ' + 
P e = V, which once projected on each component K is 
just the complementarity relation Eq. (j 16[) . 



Appendix C: Detailed balance relations 

This appendix is dedicated to proving the detailed bal- 
ance relations for IZi^j and Bi. Defining the contribu- 
tions of individual states to the partition function, 



Q K ee g K e- EK ' T * 



(CI) 



and similarly for Q.;, we make use of the standard prin- 
ciple of detailed balance for rates connecting individual 
states, 



{K 



Rk^>l — QlRl^k, 



(C2) 



and similarly QkRk^i — QiRi^K- 

We begin by defining the N x N nonsingular diagonal 
matrix F that is proportional to the equilibrium abun- 
dances, 



(C3) 



Frl = QkSkl- 



Then Eq. (|C2|) combined with the definition Eq. (|23|) 
implies that FM is symmetric. It therefore follows that 
its matrix inverse M _1 F _1 is symmetric, and hence that 



(M" 1 )^ (M- 1 )^ 



Ql 



Qk 



(C4) 



The transition rate, Eq. (fT^jl . can be expanded using 
Eq. (HH) as 

Ki^j = Ri^j + YsCM-^KLRi^KRL^j. (C5) 

K,L 



We then see that: 



K,L 



QiRi^j + Qjf(M 1 )klRk^iRl 

K,L 

QjRj^i + $3 QUM-^lkRk^Rl 

K,L 

QjRj^i + Yj 1 )LKR]<^iRj^L 



-+3 



K.L 



(C6) 



where we have used Eq. (|C2[) twice and in the third equal- 
ity used Eq. ([U4]), This proves Eq. (l2"Tj) . 
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We may also relate the effective recombination and 
photoionization rates. To do so, we consider the case 
of T m = T r and define 



/ 27T ( U e T r X " 



(C7) 



Then Eq. ((3]) can be written as Q n i(3 n i = qctni- Using 
Eq. (JS7D, we see that 



K,L 

= qa>i + J2 QkRk^M-^klPl 

K,L 

= qa>i + QlRk^M-^lkPl 

K,L 

= qai + ^ qRK^i(M~ 1 ) LK a L 

K,L 

= qA, (C8) 

where in the last equality we have used Eq. (|17j) with the 
P l K determined by Eq. <j26j). This proves Eq. (1221) . 



history using perturbation theory around the Saha ap- 
proximation, which we describe here. The idea is that the 
actual ionization fraction x e {z) is slightly greater than 
the Saha equation would predict because there are more 
recombinations than ionizations (i e < 0) and a slight 
deviation from thermodynamic equilibrium is required 
in order to drive this imbalance. We may thus take an 
ODE for the recombination history (for simplicity we use 
the Peebles ODE with an updated recombination co- 
efficient ae [iH), and Taylor-expand it around the Saha 
ionization fraction: 



i£(x e ,z) = D x {x e - xT ha ) + 0(x. 



„Saha\2 



y. (di) 



Here the superscript denotes the Peebles ODE, 

and the zeroeth-order coefficient in the Taylor series 
vanishes since thermal equilibrium considerations imply 
ig (xf aha , z) = 0. The coefficient D\ may be obtained by 
numerical differentiation; we use a two-sided finite differ- 
ence with Ax e = ±0.01xf aha (l - xf aha ). Then we obtain 
a post-Saha corrected solution by setting the left-hand 
side of Eq. (|Dlj) with if 3, ha (z) (again obtained by a two- 
sided finite difference with Az — ±1): 



Appendix D: Implementation of post-Saha 
correction at early times 



^Saha 



•Saha 



(D2) 



At the highest redshifts, the ODE describing hydrogen 
recombination is stiff, and we follow the recombination 



At the transition redshift z t = 1570, we use x c °"(z t ) as 
an initial condition for the integration with our full ODE. 
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